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Abstract 

The problem of image reconstruction in thermoacoustic tomogra- 
phy requires inversion of a generahzed Radon transform, which inte- 
grates the unknown function over circles in 2D or spheres in 3D. The 
paper investigates implementation of the recently discovered backpro- 
jection type inversion formulas for the case of spherical acquisition in 
3D. A numerical simulation of the data acquisition with subsequent 
reconstructions are made for the Defrise phantom as well as for some 
other phantoms. Both full and partial scan situations are considered. 
The results are compared with the implementation of a previously 
used approximate inversion formula. 

1 Introduction 

The idea of thermoacoustic tomography (TCT, sometimes also called TAT) ^3 
[HI Cni EH EHj can be briefly described as follows (see Figure [l]). A short pulse 
of radiofrequency (RF) electromagnetic waves is sent through a biological ob- 
ject heating up the tissue. It is known that the cancerous cells absorb several 
times more RF energy than the healthy ones P^. As a result a significant 
increase of temperature occurs at the tumor locations causing a thermal ex- 
pansion of cancerous masses pressing on the neighboring healthy tissue. The 
created pressure wave is registered by the transducers located on the edge of 
the object. Assuming speed of propagation of these acoustic waves constant 
inside the object (an assumption which is not always correct, but is satis- 
factory, e.g. for mammography), the signals registered at any transducer 
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Figure 1: Sketch of a TCT system 



location are generated by the inclusions laying on a sphere centered at that 
location. In fact the measured data are the integrals of RF absorption coeffi- 
cient / over those spheres or, in other words, the spherical Radon transform 
Rf of the RF absorption coefficient /. Hence, to reconstruct the image one 
needs to invert the spherical Radon transform. 

Definition 1.1. The spherical Radon transform of f is defined as 



where da{x) is the surface area on the sphere \x — p\ = r centered at p & M.^. 

In the definition we allow arbitrary set of centers p and radii r. However, 
from the dimensional consideration it is clear, that this mapping is overde- 
termined. The tomographic motivation suggests to restrict the set of centers 
to a surface S" C (the set of transducers' locations), while not imposing 
any restrictions on the radii. In this paper we will deal only with the case of 
spherical acquisition (i.e. the transducers are located on a unit sphere) and 
from now on we will suppose \p\ = 1. 

Two different approaches have been used to derive exact inversion for- 
mulae for this case. Fourier- Bessel and spherical harmonic expansions result 
in solutions written as an infinite series for two and three dimensions re- 
spectively ^HIC2]- The TCT analog of p-filtered backprojection inversion is 
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derived in P 

as well filtered backprojection type version 

Both formulas can be generalized to higher odd dimensions [0]. Notice that, 
as one can expect for a codimension 1 Radon transform in 3D, the formulas 
are local. 

In Section 2 we describe the numerical simulation of the data acquisi- 
tion. The reconstruction algorithms based on the p-filtered backprojection 
formula ((H) and the filtered backprojection one © are discussed in Section 3. 



2 Data Simulation 

The region of reconstruction is the unit ball centered at the origin (see fig. 121) ■ 
All phantoms considered in the paper are sums of indicator functions sup- 
ported in ellipsoids completely contained inside the unit ball. The transduc- 
ers are located on the surface of the unit sphere. 




Figure 2: Phantom setup 



We parameterize the transducer location by two angles {(pTidx)-, where 
(j)T G [0, 27r) is the the azimuthal angle in the x7/-plane and 6^ G [0, vr] is the 
polar angle measured from the z-axis. 
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Since the spherical Radon transform is hnear, it is enough to create pro- 
jections for phantoms with a single ellipsoid and then superimpose the projec- 
tions. For a single ellipsoid the data measured at a fixed transducer location 
at a given moment (i.e. for fixed {(j)T,0T,'r)) is the surface area of a part of 
the sphere of integration cut by the intersecting ellipsoid. It can be expressed 
ClS db finite sum with terms of the form 
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(3) 



where each such term corresponds to a connected component of the intersec- 
tion. Here and 6 parameterize the sphere of integration and are indepen- 
dent of (pT and 6t, which parameterize the transducer location. The angles 
6i{4>) and ^2(0) are defined by the intersection of the integration sphere and 
the phantom's ellipsoid. The cosines of these angles can be found from the 
solution of a quartic equation describing that intersection. 

In the numerical results presented below the quartic equation is solved 
using the MATLAB built-in function "roots" . By adding up these roots in an 
appropriate way we obtain the inner integral with respect to the polar angle 
6 in equation Q. The result is a function of azimuthal angle 0, which we 
will denote F{(j)). Depending on the location and parameters of the ellipsoid, 
F{(j)) might be either a smooth vr-periodic function of 0, or a piecewise smooth 
one (see fig. EI)- In the first case we compute its values at uniformly discretized 
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Figure 3: F{(f)) for two different locations of an ellipsoid 

locations on the interval [0, vr] and use the trapezoidal rule to compute the 
integral. For F(0) € C^, numerical integration using the trapezoidal rule is 
accurate to 0(Ji^) If, however, -F(0) is only piecewise smooth on [0,7r], 
then we locate the pieces of supp F{(f)) where it is smooth and use Gaussian 
quadrature to integrate over each piece. 
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3 Reconstruction 



Once we have generated the projection data, we reconstruct the original 
indicator functions of the phantoms. The reconstruction algorithms are based 
on the p-filtered backprojection ((H) or the filtered backprojection ((21) . 

The integrals over the unit sphere in (0) and ((2)) are computed as double 
integrals with respect to the azimuthal angle (pT and the polar angle 6t- The 
function to be integrated is periodic with respect to 0t making the trape- 
zoidal rule an appealing quadrature choice. Integration with respect to 9t is 
done by Gaussian quadrature. The Laplace operator is implemented through 
the Matlab built-in function "del2" . The reconstructions were generated us- 
ing Matlab 5.0. 

In the results below the resolution is 256 x 256 x 256 over a 2 x 2 x 2 
volume resulting in isotropic pixel dimension of 1/128. 

The algorithm is tested on the Defrise phantom which consists of five thin 
ellipsoids symmetrically centered along the 2;-axis (see fig. |3)). We numerate 
them from 1 to 5 starting with the lowest. 



ellipse number 


center = (xq, 2/o, 2^0) 


semiaxes lengths = (ea;,ey,e^) 


1 


(0,0,-0.64) 


(0.65,0.65,0.08) 


2 


(0,0,-0.32) 


(0.85,0.85,0.08) 


3 


(0,0,0) 


(0.9,0.9,0.08) 


4 


(0,0,0.32) 


(0.85,0.85,0.08) 


5 


(0,0,0.64) 


(0.65,0.65,0.08) 




Figure 4: The Defrise phantom slice along the plane y=0 
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3.1 Full scan data 



The data was acquired from the transducers located discretely over the sphere 
in the following way. The azimuthal angles of the transducer locations were 
uniformly discretized to iV^ = 400 points between and 27r. The polar angles 
of the transducer locations corresponded to Nq = 200 Gaussian nodes on the 
interval form to vr, as described in the previous section. The radii of the 
integration spheres were uniformly discretized to Nr = 200 points from to 
2. The reconstruction was done by both methods: filtered backprojection 
(FBP) and p-filtered backprojection. 
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(a) FBP 



(b) p-filtered reconstruc- 
tion. 
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Figure 5: Reconstructions and profiles of the Defrise phantom along the 
center x = slice. Dashed lines correspond to the center x = = y profile; 
solid lines correspond to x = 0,y = OA 
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The obtained results validate reconstruction formulas (P) and (j2I) (see 
fig- El)- In both cases the Defrise phantom has a good reconstruction every- 
where except along the z-axis {x = y = 0), where some noise is present which, 
while not always noticeable on reconstructions, is visible on the graphs. The 
reason for appearance of that noise is the correlation of numerical errors 
along that axis of phantom's symmetry and is discussed in section EH 



3.2 Partial scan data 

Half-scan reconstructions were done using data from only the eastern hemi- 
sphere {N^ = 200, Ng = 200) or the southern hemisphere {N^ = 400, 
Nq = 100). These hemispheres are highlighted in Figure IHl The rest of 
the data has been zero-filled. 




(a) p-filtered BP; 1/2- 
scan with respect to (fi. 
The profiles along the 
lines are plotted below. 



(b) FBP; 1/2-scan with 
respect to 9. The pro- 
file along the solid line 
is plotted below. 
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Figure 6: Partial scan reconstructions of the Defrise phantom. 
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It is known j2ni UHl UH 1201 1221 120] that in case of incomplete data one can 
expect to recover stably only certain parts of the image the rest of it being 
blurred out. Namely some parts of the wavefront set of the image will be lost. 
For our phantom the singularities are jump discontinuities (edges) of imaged 
value / across an interface / (a surface of an ellipsoid in 3D). The wavefront 
WF{f) of / in this situation is the set of pairs (x,n), where a; is a point 
on J, and n is a vector normal to J at x. As it was shown in ^3 121] using 
microlocal analysis, a point {x,n) G WF{f ) can be stably detected from the 
Radon data, if and only if Rf includes data obtained from a sphere passing 
through X and normal to n. In other words, one can see only those parts 
of an interface, that can be tangentially touched by spheres of integration 
centered at available transducer locations. The rest of the interface will be 
blurred. 

Edges in the Defrise phantom were reconstructed in Figure as expected. 
When the data is collected from the eastern hemisphere there are enough 
spheres to touch tangentially all edges in the eastern hemisphere (see fig. IIJ) 
but none to do it in the western hemisphere. That is why the locations of the 
edges in the eastern hemisphere were correctly reconstructed while those in 
the western part were blurred. When the data is collected from the southern 
hemisphere there are enough spheres to touch tangentially all edges in the 
Defrise phantom, hence all of them were resolved. 

From the geometric description above it is not hard to see that there may 
exist certain regions of reconstruction (locations of x, sometimes called audi- 
ble zones) where any possible pair (x, n) belonging to WF{f) is recognizable 
from Rf. In our examples, when the data is collected from the eastern or 
southern hemisphere, these regions are the eastern and southern half of the 
unit ball correspondingly. 

Notice that the image values were not reconstructed correctly, since part 
of the data was missing. However, certain iterative techniques allow one to 
improve substantially the image values in the audible zone [221 1211 121] • 

3.3 Comparison with an approximate backprojection 

In early experimental work on thermoacoustic tomography, an approximate 
backprojection formula was used. It was written in analogy with the back- 
projection of regular Radon transform and looked similar to equation (^, 
except the missing weight factor ^^z^- The composition of this operator 
with the direct Radon transform is an elliptic pseudo-differential operator 
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of order zero (see, e.g., [IZl HUl HHI-) Thus the locations and "strengths" of 
image singularities are recovered correctly. However the values of the image 
function will not be recovered correctly. The obtained reconstructions val- 
idate the predictions correctly recovering locations of edges. The values of 
image functions are accurate near the center where r ~ 1 but degrade slowly 
with distance from the origin, as expected. 




(a) Reconstruction of 
Defrise phantom at 
256 X 256 resolution 
along the center x = 
slice 



(b) Sphere of radius 
0.7 recon'd at low- 
res 64 X 64. Be- 
low: center profiles 
y — on three different 
slices. Solid x = ±0.1, 
dashed x = ±0.3, 

dash-dot x ~ ±0.5 





Figure 7: Approximate FBP shows low-frequency shading. 



9 



3.4 Errors in reconstruction 



As it was mentioned before, the reconstructions of Defrise phantom have 
some noise along the axis of phantom's symmetry x = y = (see figs. El)- 
To discuss the reasons of appearance of that noise we consider reconstructions 
of some simpler phantoms consisting of indicator functions of a perfect ball. 
This allows us to compute the Radon transform analytically, hence to exclude 
the errors in the data simulation. For every fixed po, Rf{po,r) is a third 
order polynomial with respect to r for 0<ri<r<r2<l and is zero for 
every other r. Filtered backprojection requires differentiating with respect 
to the radial variable r. We used centered finite differences to estimate 
d? / dr'^Rf{p,r) which is exact on the third degree polynomials. Therefore, 
we compute (P / dr'^Rf{p^r) exactly for all radii, r, at least 2Ar away from 
ri and r2. Hence the only errors in numerical differentiation that spread 
into the backprojection come with the data from spheres close to the ones 
touching tangentially the phantom ball. None of these spheres passes inside 
the phantom ball, hence backprojection at those points is free of errors from 
numerical differentiation (see fig. |H1). 
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Figure 8: FBP errors along the axis of symmetry. 

Now let us consider a point pi on the axis of symmetry of the ball phan- 
toms (the line connecting the center of the phantom ball and the origin). 
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There are two sets of spheres that pass through that point and touch the 
phantom ball tangentially. The spheres in the first set contain the phantom 
ball, while the spheres in the second set do not. A 2D slice of this scenario 
is presented in Figure M 



Figure 9: The spheres C3 and C4 contain the phantom ball, while CI and 
C2 do not 

Notice that all spheres in the same set have the same radius. So the errors 
from the numerical differentiation that they will bring into the backprojection 
algorithm are absolutely the same. The axis of symmetry is the only location 
in the reconstruction region where these errors are perfectly correlated. This 
resonance increases the magnitude of errors resulting in the noise along the 
symmetry axis on reconstructed images (see fig. |H1)- 

In case of ellipsoids in the Defrise phantom everything said above holds. 
In fact magnitude of errors is five times bigger since there are five ellipsoids 
with the same axis of symmetry there. At the same time, the reconstruction 
of an ellipsoidal phantom without any rotational symmetry has no axis of 
emphasized errors (see fig. ITUj). 
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50 100 150 200 250 50 100 150 200 250 



Figure 10: An ellipsoidal phantom with center at (0,0.2,-0.1) and semiaxes 
lengths equal to (0.4,0.3,0.5). 



4 Conclusion 

We have implemented a straightforward numerical validation of both FBP 
and p-filtered inversion formulae for TCT data on the high-frequency Defrise 
phantom. FBP and p-filtered have virtually identical performance with noise- 
free simulated data from this high-contrast object. Artifacts due to numerical 
errors are more severe in FBP than p-filtered images and might be reduced by 
mollification techniques ^T]. Comparing FBP and p-filtered performance in 
the presence of noise and for low-contrast detectability will appear in future 
pubhcations. 
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